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Abstract -We study the heterogeneous nucleation of Ising model on complex networks under 
a non-equilibrium situation where the impurities perform degree-biased motion controlled by a 
parameter a. Through the forward flux sampling and detailed analysis on the nucleating clusters, 
we find that the nucleation rate shows a nonmonotonic dependence on a for small number of 
impurities, in which a maximal nucleation rate occurs at a = 0 corresponding to the degree- 
uncorrelated random motion. Furthermore, we demonstrate the distinct features of the nucleating 
clusters along the pathway for different preference of impurities motion, which may be used to 
understand the resonance-like dependence of nucleation rate on the motion bias of impurities. 
Our theoretical analysis shows that the nonequilibrium diffusion of impurities can always induce a 
positive energy flux that can facilitate the barrier-crossing nucleation process. The nonmonotonic 
feature of the average value of the energy flux with a may be the origin of our simulation results. 


Introduction. — Nucleation is an activated process 
which initiates the decay of a metastable state into a more 
stable one jl] driven by fluctuation. Many important dy¬ 
namical processes on real-world scenarios, such as crys¬ 
tallization 0[3], fractures HE], glass formation [|6j, and 
protein folding [7], to list just a few, are concerned with 
nucleation. For many decades, our understanding of nu¬ 
cleation has been dominated by the classical nucleation 
theory (CNT), and it has been applied not only to the 
liquid-gas and liquid-solid systems, but also to regular lat¬ 
tices in Euclidean space. For instance, in two-dimensional 
lattices, Allen et al discovered that shear can enhance 
the nucleation rate and the rate peaks at an intermedi¬ 
ate shear rate [8] . Sear found that a single impurity may 
considerably enhance the nucleation rate [5]. Page and 
Sear reported that the existence of a pore may lead to 
two-stage nucleation, and the overall nucleation rate can 
reach a maximum level at an intermediate pore size ma¬ 
in three-dimensional lattices, the nucleation pathway of 
the Ising model has also been studied by Sear and Pan 
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mm- In addition, the validity of CNT has been tested 
in other Euclidean space mm- 

Since many real systems can be properly modeled by 
network-organized structure [20H22], it is thus an inter¬ 
esting topic to explore nucleation process in complex net¬ 
works. Recently, we have studied nucleation dynamics of 
the Ising model on scale-free (SF) networks [23], Erdos- 
Renyi (ER) networks [23] and modular networks [25]. We 
found that, for homogeneous nucleation on SF networks, 
many small isolated nucleating clusters emerge at the early 
stage of the nucleation process, until suddenly they form 
the critical nucleus through a sharp merging process, and 
the nucleation rate decays exponentially with network size. 
For homogeneous nucleation on ER networks, there always 
exists a dominant nucleating cluster to which relatively 
small clusters are attached gradually to form the critical 
nucleus. For modular networks, as the network modular¬ 
ity worsens the nucleation undergoes a transition from a 
two-step to one-step process and the nucleation rate shows 
a nonmonotonic dependence on the modularity. For het¬ 
erogeneous nucleation, target impurities are shown to be 
much more efficient to enhance the nucleation rate than 


P-1 





C. Shen et al. 


random ones. However, in our previous work, impurities 
are considered to be fixed in some nodes with the new 
phase. As we know, mobility is a ubiquitous feature of real 
systems [26H30j ■ including the mobility of impurities, and 
may drastically influence the dynamical evolution. For ex¬ 
ample, it has been reported that mobility promotes syn¬ 
chronization mm, enhances signal response [32] , affects 
contagion processes [34], and tunes biodiversity [35] . In 
addition, impurities may be caused by the vacancy de¬ 
fects with no interaction, not with new phase. How the 
impurity motion would influence the nucleation rate and 
pathway is still an open question. Motivated by this, we 
will study the different roles of the motion in the forma¬ 
tion of nucleating clusters, which can reveal the nucleation 
pathways of the Ising model in the underlying networks. 

In the present work, we adopt the recently proposed for¬ 
ward flux sampling (FFS) [36] approach, which is efficient 
and easy to implement to study rare events, and employ 
SF networked Ising model. Ising model is a paradigm for 
many phenomena in statistical physics and widely used 
to study the nucleation process. By introducing degree- 
biased random walks for impurities on the network, we 
find that for small number of impurities the nucleation 
rate shows a nonmonotonic dependence on the bias pa¬ 
rameter of the motion of impurities, in which a maximal 
nucleation rate occurs at the situation where the impu¬ 
rities perform random motions. Furthermore, we show 
that there are different properties of the nucleating clus¬ 
ters along the pathway corresponding to different impuri¬ 
ties bias motions. 

Model and method. — 

Model. We consider the non-equilibrium Ising model 
with mobile impurities on complex networks consisting of 
N normal nodes and w defect nodes called impurities. 
Each normal node is endowed with a spin variable Sj that 
can be +1 (up), or —1 (down), and each defect node is 
endowed with spin 0 (impurity). The Hamiltonian of the 
system is given by 

H — J ^ ] A ij SiSj h ^ ] Si , (1) 

i<j i 

where J is the coupling constant and h is the external 
magnetic field. For convenience, we set J = 1 in the fol¬ 
lowing discussions. The elements of the adjacency matrix 
of the network take Ajj = 1 if nodes i and j are connected 
and Aij = 0 otherwise. The degree, that is the number of 
neighboring nodes, of node i is defined as ki = ]Cj!=i Ajj. 
Notice that, without defect nodes, there exist a number of 
simulations and analytical results for the Ising model in 
ER and SF networks EZHU]- 

The dynamical evolution of our model has two ingredi¬ 
ents: Spin-flip and Impurity diffusion. In each time step, 
we attempt to perform the following two types of dynam¬ 
ics. 1) Spin-flip: we randomly chose a normal node and 
attempt to flip its spin according to the Metropolis accep¬ 
tance probability min(l, e -/3AE ) [42], where /3 = 1 /(fc^T) 


with the Boltzmann constant fcs and the temperature T, 
and A E is the energy change due to the flipping process; 2) 
Impurity diffusion: After a spin of a normal node has been 
attempted to be flipped, we then randomly choose a de¬ 
fect node(impurity node) i and exchange the spin Si with 
that Sj of one of the nearest neighboring normal nodes j 
according to the probability 




D 


E 




( 2 ) 


Here D is the diffusion constant, the sum is taken over 
all the nearest neighboring normal nodes of i, and a is 
a tunable parameter which biases the impurities’ motion 
either towards low-degree nodes (a < 0) or towards hubs 
(a > 0). For a = 0, we recover the standard (unbiased) 
random walk. 

In general, with the increment of T the system under¬ 
goes a second-order phase transition from an ordered state 
to a disordered one at the critical temperature T c . Below 
T c the system prefers to be in a state with all spins up 
or down. Given a small external field, one of these two 
states will become metastable, and if initiated predomi¬ 
nantly in this metastable state, the system will remain for 
a significantly long time before it undergoes a nucleation 
transition to the thermodynamically stable state. We are 
interested in the rate and pathways for this transition. 

FFS method. The FFS method has been successfully 
used to calculate rate constants and transition paths for 
rare events in equilibrium and non-equilibrium systems 
[SHHBIMlElil]- This method uses a series of interfaces 
in phase space between the initial and final states to force 
the system from the initial state A to the final state B in 
a ratchetlike manner. First, we define an order parameter 
A(:r), where x represents the phase-space coordinates, such 
that the system is in state A if A(;r) < Ao and state B if 
A (a;) > Am, while a series of nonintersecting interfaces 
Ai(0 < i < M) lie between states A and B, such that 
any path from A to B must cross each interface without 
reaching A^+i before A j. The transition rate R from A to 
B is calculated as 

n M-l 

P(A i+1 |Ai) (3) 

z=0 

where i>A,o is the average flux of trajectories crossing Ao 
in the direction to B. P(Am|Ao) = Tli^o 1 P (A»_(_i|A*) is 
the probability that a trajectory crossing Ao in the direc¬ 
tion to B will eventually reach B before returning to A, 
and P (Ai+i|Aj) is the probability that a trajectory which 
reaches A i, having come from A , will reach A^+i before 
returning to A. For more detailed descriptions of FFS 
method, please see Ref. [45] . 

In this work, we perform Monte Carlo simulation and 
use FFS to study nucleation rate and pathways of the non¬ 
equilibrium phase from the metastable spin phase. Specifi¬ 
cally, we set T < T c , h = 0.5 and start from an initial state 
with s = —1 for most of the spins. We define the order pa¬ 
rameter A as the total number of up spins in the network. 
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Fig. 1: (Color online) The nucleation rate Ini? as a function of 
a for different values of the diffusion rate D. The dotted line 
indicates the result without impurities, i.e., w = 0. Parameters 
are N = 1000, the average network degree (fe) = 6, w = 2, 
h = 0.7, T = 2.59, A 0 = 130 and Am = 880. 

The spacing between adjacent interfaces is fixed at 3 up 
spins. We perform 1000 trials per interface for each FFS 
sampling, from which at least 200 configurations are saved 
in order to investigate the statistical properties along the 
nucleation pathway. The results are obtained by averag¬ 
ing over 20 independent FFS samplings and 50 different 
network realizations. 

Results and Discussion. — 

Nucleation rate. In what follows, we employ a 
Barabasi-Albert (BA) SF network, whose degree distri¬ 
bution follows a power law p{k) ~ fc -7 with the scaling 
exponent 7 = 3 [46! . 

To begin, we fix the small number w = 2 of impurities 
and vary the diffusion constant D to investigate how the 
nucleation rate R (in unit of MCstep -1 spin -1 ) evolves 
with controlling parameter a. Figure |T] shows the depen¬ 
dence of the logarithm of the nucleation rate In R on a for 
different values of D. One can observe that In R exhibits a 
resonance-like behavior with the increment of a. That is, 
there exists an optimal value of a at a = a op t, correspond¬ 
ing to the maximum R. Interestingly, this phenomenon is 
robust against the diffusion constant D. This result in¬ 
dicates that random motions of impurities corresponding 
to a op t = 0 is more favorable to nucleation than degree- 
biased motions. In addition, for any given values of a we 
find that In R increases monotonously with D , indicating 
that impurities mobility is always in favor of nucleation. 
That is, the larger the mobility rate D is, the larger In R 
becomes. The dotted line indicates the result without im¬ 
purities, i.e., w = 0. Obviously, for any given value of D 
impurities may considerably enhance the nucleation rate, 
which is consistent with [5]. 

It is worthy noting that this nontrivial dependence is 


unobservable if the number of impurities w becomes rela¬ 
tively large. Figure [2] shows the dependence of In I? on 
a for different values of w. Clearly, for small number 
of impurities, say ui = 1,3, 5, 8,10, one can always ob¬ 
serve an interesting mobility induced resonance-like be¬ 
havior in accordance with Fig.l. While for big w , say 
w = 15, 20, 25, 30, In R increases monotonously with a for 
a < 1 and then approaches a constant value for a > 1. 
Other values of w have also been investigated; the quali¬ 
tative results are the same and not shown here. But for 
w = 0 indicated by the dotted line, i.e.,without impurities, 
In I? is considerably less than that of impurities. 



Fig. 2: (Color online) In A as a function of a for different values 
of w. The dotted line indicates the result without impurities, 
i.e., w = 0. Other parameters are same as in Fig.l except for 
D = 0.5 and h = 0.5. 

Nucleation pathway. To elucidate the detailed char¬ 
acteristics along the nucleation pathway for different bias 
motions, we save lots of configurations generated by FFS 
and perform detailed analysis on the nucleating clusters, 
including the relative size of the largest and the second 
largest cluster, average degree of the cluster nodes and 
the number of clusters. According to CNT, there exists 
a critical nucleus size A c of the new phase, above which 
the system grows rapidly to the new phase. Herein, we 
mainly focus on the nucleation stage where A < A c . In our 
simulation, we determine A c by computation of the com¬ 
mittor probability Pb , which is the probability of reach¬ 
ing the thermodynamic stable state before returning to 
the metastable state. As commonly reported in the litera¬ 
ture mm, the critical nucleus appears at Pb(A c ) = 0.5. 
Since A c are different for different bias parameters, we thus 
introduce A/A c as the control parameter. 

Following the previous study [24], we introduce the rel¬ 
ative size S m axi Ssec of the largest and the second largest 
nucleating cluster, which are defined as the ratios of the 
number of up spins within the largest and the second 
largest cluster to the total number of up spins respectively. 
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Fig. 3: (Color online) (a) The relative size Smax, S sec of the 
largest and the second largest nucleating cluster respectively, 
as a function of A/A c . (b) The average degrees kmax, k se c of the 
nodes within the largest and the second largest nucleating clus¬ 
ter respectively, as a function of A/A c . (c) and (d) correspond 
to the average degree (/:„„) of the nodes inside nucleating clus¬ 
ters and the number n c of nucleating clusters respectively, as a 
function of A/A c for a = —3, 0, 3, corresponding to motion pre¬ 
ferring to low-degree nodes, to random nodes and high-degree 
nodes respectively. Symbols for different motion bias in (b), 
(c) and (d) are same as in (a). Other parameters are same as 
in Fig.l except for D = 0.5 and h = 0.5. 


Smax and S sec (averaged over the ensemble at each inter¬ 
face) as a function of A/A c are plotted in figure [3] (a). 
Clearly, one can see that S max for a = 0 (solid red cir¬ 
cles) is always larger than those for a = —3 (solid black 
squares) and a = 3 (solid blue triangles). Specifically, 
at A/A c = 0.35, S max grows already more than 50% for 
a = 0, while it is less than 30% for a = —3 and about 
20% for a = 3, as shown by the dashed gray lines in fig¬ 
ure (Ha). But when A/A c = 1 they almost tend to 100% 
together. This difference means that for unbiased random 
motion S ma x grows fast at the very beginning following by 
a relatively slow increasing, while for biased motion, S max 
increases slowly at first and then rapidly when approach¬ 
ing the critical nucleus. From figure [3] (a) one can also 
observe that the relative size S sec (denoted by the empty 
symbols) is greatly less than 5 mal , indicating that the 
nucleation is dominated by the largest nucleating cluster. 

We also plot (k m ax) and ( k sec ), defined as the aver¬ 
age degrees of the nodes within the largest and the sec¬ 
ond largest nucleating cluster respectively, as a function 
of A/A c in figure [3] (b). Clearly, (kmax) for a = 0 indi¬ 
cated by the solid red circles, is always larger than those 
for a = —3 and a = 3 indicated by the solid black squares 
and solid blue triangles respectively. Strikingly, at the 
early nucleating stage, ( k max ) grows sharply for the for¬ 
mer, while it grows gradually for the latter. Furthermore, 



Fig. 4: (Color online) Schematic illustration of an impurity 
diffusion, (a) indicates an impurity node occupied a fc-degree 
node (filling node) will exchange its present position with a spin 
on a fc'-degree node (arrow node), and (b) denotes the reverse 
process, m' and m" denote the average magnetization of a ran¬ 
domly chosen neighboring normal node and defect (impurity) 
node respectively. About their definitions, please see Eos. (14151) . 


it is found that ( k sec ) is greatly less than ( k max ), which 
suggests again that the largest nucleating cluster domi¬ 
nates the nucleation. 

In addition, we also investigate the average degree 
(knew) of the nodes inside the new phase and the num¬ 
ber n c of nucleating clusters, and plot (k new ) and n c as a 
function of A/A c in figure [3] (c) and (d) respectively. As 
shown, (k new ) increases monotonically with A/A c for dif¬ 
ferent a = —3,0,3, which means the new phase tends to 
grow from those nodes with smaller degrees. Nevertheless, 
for different preference of impurities motion it shows the 
distinct features along the nucleation pathway. For a = 0, 
(knew) grows fast at the very beginning following by a rela¬ 
tively slow increasing. For a = —3 and 3, (k new ) increases 
slowly at first and then grows fast until approaching the 
critical nucleus. Such a scenario is consistent with figures 
[3] (a) and (b). From figure [3] (d) one can observe that 
n c shows non-monotonically dependence on A/A c for dif¬ 
ferent a. Especially, the number of clusters for a = 0 is 
always less than that for a = —3,3. On the other hand, 
n c approaches the same magnitude near the formation of 
critical nucleus for three different a, which suggests that 
it is easier for the critical nucleus comes into being for un¬ 
biased random motion than that for others. This result is 
also consistent with the picture shown in figures [3] (a) to 
(c). 

To further understand the nontrivial effect of the 
nonequilibrium diffusion on nucleation, we will evaluate 
the average energy change due to the diffusion process be¬ 
tween an impurity and a spin. To begin, let us consider 
such a process on a link connecting a fc-degree node and a 
fc'-degree node. As shown in Fig|4](a), an impurity occu¬ 
pied fc-degree node will exchange its present position with 
a spin of fc'-degree node. The reverse process is shown 
in FigQJb). Let qk and nik denote the probability of a 
fc-degree node occupied by an impurity and the average 
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magnetization of a fc-degree node occupied by a spin, re¬ 
spectively. For a spin node on degree uncorrelated net¬ 
works, the average magnetization m! of a randomly chosen 
neighboring node is written as, 

m ' = ~ 9fc) (4) 

While for an impurity node, the average magnetization of 
a randomly chosen neighboring node is written as 

m " = Y. k ~ 9fc) (5) 

where the subscript k — 1 other than k is based on the 
consideration that there is no any interaction between an 
impurity node and its neighbors. For the case of FigQJa), 
the energy change is written as 

AEi = — Jm k '[(k — 1 )m" — ( k ' — l)m'] (6) 

Analogously, for the case of Figj4|b) 

AE 2 = —Jmk[(k' — l)m" — (k — l)m'] (7) 


The energy change due to a diffusing exchange taking plac¬ 
ing on a k — k' link can be expressed as 

A E = - (q k pk^k'AEi + q k 'Pk'^kAE 2 ) (8) 

where Pk^k 1 is the diffusion rate of an impurity from k- 
degree node to fc'-degree node. According to our model, 
it can be expressed as 


Pk — >k' 



k' a 


k'P(k') ]ln 
(k) ^ 


p (k)k ,a 


(9) 


Assuming the diffusion is quasi-static process that satisfies 
the detailed balance conditions, 


QkPk^k' = Qk'Pk'^k 


( 10 ) 


The requirement can lead to the expression of q k , 


wP(k)k a+1 

qk ~Z k P(k)k a+1 


( 11 ) 


Averaging over all possible links on networks, one arrives 
at the average energy change due to a nonequilibrium dif¬ 
fusion process, 


(A E) = £ Ikk’AE (12) 

k,k' 

where l k k’ = kk'P(k)P(k') / (k) 2 is the probability that 
a randomly chosen link to connect a pair of nodes with 
degree k and k'. 

Next we will calculate the average magnetization using 
heterogeneous mean-field theory. Following Ref. [47|, one 
has 


TOfc = tanh[/3 Jkm' + (3h\ (13) 



a 

Fig. 5: (Color online) The average energy change (A E) as a 
function of a. The solid line indicates the results of the theo¬ 
retical prediction, and solid symbols that of MC. Parameters 
are same as in Fig.l except for D = 1.0 and h = 0.5. 


Substituting Eq. (fl3l) with Eqs.(|4][5]), we arrive at the self- 
consistent formulations of m' and m" that can be numer¬ 
ically calculated. 

Figure [S] shows the results of (A E) as a function of 
a , where the solid line denotes the result obtained from 
Ea. (fl2l) . and the symbols that of MC simulations. Clearly, 
the theory can reproduce qualitatively well the main char¬ 
acteristic: there exists an optimal motion bias where the 
average energy change reaches the maximum. Further¬ 
more, it is found that (A E) are always larger than zero 
for any motion bias a, which indicates that the impurity’s 
mobility can always facilitate the barrier-crossing nucle¬ 
ation process, akin to the drag effect of nonequilibrium 
thermodynamic forces conjugated to the exchange of im¬ 
purities and spins. 

Conclusions. — In summary, we have studied the het¬ 
erogeneous nucleation of a non-equilibrium Ising model 
with mobile impurities on complex networks. By intro¬ 
duced a tunable parameter a, the impurities can perform 
three different bias motions: a > 0 means that the im¬ 
purities prefer to visit the high-degree nodes, a < 0 the 
low-degree nodes, and a = 0 recovers to completely ran¬ 
dom motion. Interestingly, it is found that the nucleation 
rate is not a monotonic function of a for small num¬ 
ber of impurities, i.e., there exists an optimal value of 
a = 0, leading to the fastest nucleation rate. Especially, 
the optimal value of the controlled parameter does not 
change with the variation of the diffusion rate. To qual¬ 
itatively understand the underlying mechanism of such a 
phenomenon, we have performed heterogeneous mean-held 
analysis. Furthermore, we have used the FFS method and 
analyzed the nucleating clusters, and found that for differ¬ 
ent preferences of impurities motion, the nucleating clus¬ 
ters show the distinct features along the pathways. On 
the one hand, for random motion, the largest nucleating 
cluster dominates the nucleation, and the average degree 
of the nodes inside nucleating clusters grows rapidly at 
first, while for motion to the high-degree nodes or to the 
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low-degree nodes, they grow slowly at the very beginning 
following by a relatively fast increasing. On the other 
hand, the number of nucleating clusters for the former is 
less than that of the latter, especially they decreases to the 
same magnitude at the formation of the critical nucleus. 
These distinct features may mean different microscopic 
mechanisms driving the system towards nucleation. Since 
heterogeneous nucleation is essential for many dynamical 
processes on real-world scenarios, and mobility is a ubiq¬ 
uitous feature of real systems, our study may provide a 
valuable understanding for many non-equilibrium phase 
transitions taking place in networked systems and for ef¬ 
fective controlling strategy to the rate of such processes. 
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